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Preface 


The  two  parts  of  this  report  present  an  introduction  to  the  design  and  usage  of 
the  NFEARS  system  developed  jointly  at  the  Universities  of  Maryland  and 
Pittsburgh.  The  acronym  stands  for  "Nonlinear  Finite  Element  Adaptive 
Research  Solver"  and  refers  to  an  experimental  software  system  for  the  solution 
of  a  class  of  nonlinear,  stationary  boundary  value  problems  in  two  space 
dimensions  involving  two  real  parameters.  The  overall  concept  of  the  system  is 
similar  to  that  of  the  earlier  system  FEARS  (see,  e.g.,  [1 6],[1 7],[26])1  for  linear 
elliptic  problems  and  retains  many  of  the  features  of  that  system. 

The  eight  sections  of  Part  I.  describe  the  mathematical  background  of  the 
processes  incorporated  into  NFEARS,  and  the  six  sections  of  Part  II.  serve  as 
User's  Manual  for  the  program. 

As  the  earlier  program,  NFEARS  was  developed  strictly  as  a  research  tool  and 
is  expected  to  be  subject  to  continual  modifications. 


1  All  references  are  collected  in  Section  1.9  of  Part  I. 


Introduction 


As  noted  in  the  Preface,  NFEARS  is  an  experimental  software  system  for  the 
solution  of  a  class  of  nonlinear,  stationary  boundary  value  problems  in  two 
space  dimensions  involving  several  parameters.  The  specific  form  of  the 
problem  class  that  can  be  handled  by  NFEARS  will  be  given  in  the  next 
Section.  Generically,  these  are  equations  of  the  form 


F(z,X)  =  0 


involving  a  particular  nonlinear  mapping  F:  X  -4  Y.  Here  Y  is  a  suitable  function 
space  and  X  is  the  product  X  =  Z  x  A  of  another  function  space  Z  (the  state 
space)  and  a  finite-dimensional  parameter  space  A.  The  basic  formulation  of 
the  NFEARS  problems  incoporates  a  seven-dimensional  parameter  space,  but 
for  the  calculation  only  a  two-dimensional  subspace  A  is  used.  Parameterized 
nonlinear  equations  arise  in  numerous  areas.  The  applications  we  had 
especially  in  mind  derive  from  nonlinear  structural  mechanics,  but  the  system, of 
course,  is  more  generally  applicable. 

In  general,  the  set  of  all  solutions  (z,X)  of  (1)  forms  a  differentiable  manifold  M  in 
the  space  X  with  a  dimension  equal  to  that  of  the  parameter  space  A  (see  e.g. 
[19]).  In  most  applications  interest  centers  not  so  much  on  computing  a  few 
solutions  (z,X)  for  specific  values  of  the  parameters,  but  rather  on  analysing  the 
form  and  special  features  of  the  manifold.  For  example,  in  the  mentioned 
structural  problems  we  may  wish  to  determine  those  solutions  where  the 
stability  behavior  changes. 

For  the  computation,  a  finite  element  approximation  is  applied  to  the  basic 
equations.  Since  the  parameter-vector  does  not  need  to  be  discretized,  the 
resulting  equations  then  have  the  generic  form 

Fh(ZhA)  =  0 


where  now  Fh  maps  a  space  Xh  =  Zh©A  into  a  space  Yh  where  now  both  Zh  and 
Yh  are  finite  dimensional.  Since  the  parameter  space  is  unchanged,  we  may 
expect  that  the  solutions  (zh,X)  of  the  discretized  equation  form  a  differentiable 


manifold  Mh  in  Xh  of  the  same  dimension  as  the  original  manifold  M  in  X.  More¬ 
over,  we  may  imbed  Xh  into  X  which  enables  us  to  define  the  discretization  error 
as  a  suitably  specified  distance  between  M  and  Mh  (see  [12],[13]).  This  allows, 
In  turn,  for  the  design  of  a  posteriori  estimates  of  this  discretization  error  (see 
e.g.  [6], [8], [18], [20]).  As  in  the  linear  case  (see  e.g.  [2], [7], [8], [10])  such  a 
posteriori  estimates  are  not  only  important  in  providing  some  assessment  of  the 
reliability  of  the  computed  results,  but  also  in  controlling  an  adaptive  mesh- 
refinement  procedure  with  the  aim  of  obtaining  a  solution  with  appropriately 
bounded  error  behavior  using  minimal  cost  (see  [1],  [3]). 

The  basic  procedures  for  the  computational  analysis  of  the  solution  manifold  Mh 
are  the  continuation  methods,  or  incremental  methods,  as  they  are  also  called 
in  the  engineering  literature.  When  the  dimension  of  the  parameter  space,  and 
hence  of  Mh,  exceeds  one,  these  methods  require  an  a  priori  restriction  to 
some  path  on  Mh  and  then  produce  a  sequence  of  points  along  that  path.  In 
structural  problems  the  parameters  often  define  the  load  configuration  and  the 
path  is  defined  by  fixing  the  load  points  and  load  directions  and  leaving  only  a 
load  intensity  as  scalar  parameter. 

Obviously,  it  is,  in  general,  not  easy  to  develop  a  good  picture  of  a  multi-dimen¬ 
sional  manifold  solely  from  information  along  such  paths.  This  led  recently  to 
the  development  of  methods  for  the  computation  of  simplicial  approximations  of 
open  subsets  of  Mh  (see  [21], [22]). 

The  goal  of  NFEARS  is  to  provide  a  tool  for  the  above  listed  tasks.  In  fact  our 
generic  discussion  so  far  identifies  all  the  principal  capabilities  of  the  system 
which  can  be  summarized  as  follows: 

(i)  Construct  a  finite  element  discretization  of  the  given  equations  using 
biquadratic  elements  on  a  hierarchical  class  of  meshes  defined  recursively  by 
repeated  refinement  or  de-refinement. 

(ii)  With  a  given  solution  Xho  as  starting  point  and  a  path  through  Xho  specified 
by  a  selected  combination  of  the  parameters,  apply  a  continuation  process  to 
compute  a  sequence  of  points  approximating  that  path.  The  process  follows  the 


design  of  the  PITCON  system  (see  [23] ,[24]).  It  allows  also  for  the  computation 
of  specific  target  points  and  limit  points  on  the  path. 
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(iii)  At  selected  computed  points  along  the  continuation  path,  compute  a 
posteriori  estimators  for  the  discretization  error. 

(iv)  On  the  basis  of  a  "density"  and  "intensity”  concept  for  finite  element  meshes, 
use  the  error  estimators  to  control  the  modification  of  the  current  mesh. 

(v)  At  any  of  the  computed  points,  apply  an  algorithm  for  the  computation  of  a 
simplical  approximation  of  an  open  region  of  the  manifold  Mh  surrounding  the 
particular  point. 


1.1.  Problem  Class 


The  nonlinear  problems  underlying  NFEARS  are  formulated  in  a  weak  form 
requiring  the  determination  of  a  function  U  eX  such  that 

B(U,V)  =  L(V)  ,  VVeY, 


where  X,Y  are  appropriate  function  spaces  which  shall  not  be  further  identified. 
Here  B  is  a  bivariate  form 


B(U,V)  =J|j(oi'*'y.u'uvuy)Vdxdy  + 

a 

J^.x.^u.u.iyv.dxdy 

n  x  n  y 

and  L  a  linear  functional 


(<J1fx,y,U,Ux,Uy)Vydxdy 


L(V)  =  J  G2(ct2,  x,  y)  V  dxdy  +  Jg^,  x,  y)  V  dy 
n  r 

involving  the  given  functions 

(o1,x,y,U01U11U2)  ->  0(o1,x,y,U°,U1,U2)  ,  Vc^  e  R2,  x,y,  llV.U2  e  R1 
(02,x,y)  -►  G2(a2,x,y),  V  o2  e  R2,  x,y  e  R1 
(a3,x,y)  G^  (<73,x,y),  V  a3  e  R2  ,  x,y  e  R1 

In  this  formulation  we  use  the  notation: 

U  the  trial  function  in  X, 

V  the  test  function  in  Y, 

Q  a  given  domain  in  R2  which  will  be  described  in  more  detail  in 

Section  1.2, 

f  a  given  subset  of  dQ  which  will  also  be  discussed  in  Section  1.2 


x,y  the  global  coordinates  in  R2 

Y  the  arclength  of  r, 

<?i  i=1 ,2,3,  three  two-dimensional  parameter  vectors  which  will  be 

discussed  below. 

The  trial  and  test  spaces  X  and  V  are  assumed  to  allow  us  to  handle  prescribed 
Dirichlet  boundary  conditions  of  the  form 


■ ,  _  J  b(x,y) , 
u  -  \0,b(x,y), 


on  0 
on  0 


0  u  0  =  © 
1  2 


where  0  is  a  given  one-dimensional  subset  of  £2,  04  a  scalar  parameter,  and  b 
a  function  which  is  locally,  piecewise  quadratic  (see  Section  1.3). 

All  given  functions  are  supposed  to  be  sufficiently  smooth  on  their  respective 
domains.  For  the  computation  it  is  required  that  subroutines  are  available  for  the 
computation  of  the  following  functions 


90 

auk 


d20 


doj  U 


32<D 


au k  au J 


k,j  =  0, 1,2 


3G,  3G, 

G,  ,  -  ,  G.  ,  - 

2  3o2  1  dc3 

and  it  is  assumed  that  appropriate  values  are  provided  which  define  the 
boundary  function  b(x,y). 


Altogether  the  problem  involves  a  seven-dimensional  space  of  parameters 
(01,02.03.04)  e  R2  x  R2  x  R2  x  R1.  Thus  in  general,  the  set  of  all  solutions  U  is 
expected  to  be  a  seven  dimensional  differentiable  manifold  in  the  space  X  (see 
e,g,  [19]).  In  applications,  interest  rarely  centers  on  the  determination  of  a  few 
specific  solutions  of  this  problem  for  fixed  parameter-values,  but  rather  on  an 
assessment  of  the  behavior  of  these  solutions  when  the  parameters  vary;  that  is, 
on  an  analysis  of  the  form  and  features  of  the  solution  manifold.  NFEARS  is 
intended  to  provide  a  tool  for  such  an  analysis. 


Evidently,  serious  data  handling  problems  are  likely  to  arise  in  an  analysis  of  a 
seven  dimensional  manifold,  even  though,  in  principle,  the  computational 
methods,  used  here  for  approximating  segments  of  such  manifolds,  remain 
applicable.  Accordingly,  for  the  computation,  NFEARS  is  restricted  to  two- 
dimensional  sub-manifolds.  These  sub-manifolds  are  defined  by  linear 
combinations  of  the  three  parameter  vectors  <j\  ,i=1 ,2,3,  and  of  the  scalar 
parameter  <74  in  dependence  of  two  effective  parameters  X-\  and  X.2-  More 
specifically,  the  following  linear  relations  between  the  cj  and  A.j  are  assumed  to 
be  given: 

i]  Til  f0i  n  I  r,  "I 

ai  ai  Pi  0 

a  =  =  +  ,  i=1 ,2,3 

2  2  n  n2  y 

a  a.  0  p  k2 


ct4  =  a4  +  (i4 

Accordingly,  we  are  now  interested  in  the  two  dimensional  manifold  of  solutions 
(U,^i,^2)  €  X  x  R2  of  our  problem,  and  NFEARS  computes  finite  element 
approximations  of  points  along  paths  as  well  of  nodes  of  triangulations  of  open 
subsets  of  this  manifold. 


1.2.  Domains  and  Mappings 


The  domain  W  of  the  problem  is  assumed  to  be  sub-divided  by  the  user  into  a 
collection  of  sub-domains.  The  interior  of  these  sub-domains  are  called  2-D 
domains  and  will  be  denoted  here  by 

Q*  ,  k=  1 . N2 

The  closure  of  each  2-D  domain  is  a  generalized  quadrilateral  with  four  corner 
points  and  four  sides  each  of  which  can  be  either  a  straight  line  or  a  circle.  The 
(relatively)  open  line  segments  forming  the  sides  are  called  1-D  domains  (or 
lines,  for  short)  and  are  denoted  by 

aj  ,  k  =  1 

Finally,  we  call  the  corner  points  1-D  domains  (or  points  for  short)  and  denote 
them  by 

fi”,  k-l...,N0 

Thus  altogether  the  original  domain  is  decomposed  as  follows 

£2  =  (U  ft  Oj  )  u  (Uft  £2)  )  u  (Uft  £\°  ) 

The  closure  of  each  2-D  domain  is  the  homeomorphic  image  of  the  unit  square 
S  =  {(S.ti);  0  <  ^,r\  <  1  }  in  the  plane  where  all  computations  are  performed.  The 

specific  form  of  the  mapping  is  given  below.  Clearly,  these  transformations  onto 

2 

S  require  that  degeneracies  in  the  definition  of  any  2-D  domain  must  be 

2 

avoided.  In  particular,  the  angles  formed  at  the  corners  of  should  not  be  too 

close  to  0  or  k  ,  no  overlapping  sides  should  occur,  and  the  overall  shape 
should  not  be  approximately  triangular.  Figure  1.2.1  shows  some  illegal  and 
legal  configurations. 


Figure  1.2.1 


Cracks  may  be  introduced  by  using  two  1-D  domains  between  two  2-D 
domains.  If  the  crack  originates  at  the  external  boundary  of  i},  then  two  0-D 
domains  with  the  same  coordinates  should  be  used  : 


Internal  crack 


External  crack 


Q  0-D  subdomains 


Figure  1.2.2 


As  noted  before  ,  the  allowable  1-D  domains  are  either  straight  line  segments  or 
circular  arc  segments  with  a  given  radius.  More,  specifically,  any  such  domain 
0,1  is  defined  as  a  directed  arc  from  a  starting  point  with  coordinates  (x^y,)  to  an 

endpoint  with  coordinates  (x2,y2).  Of  course,  both  of  these  points  are  0-D 
domains.  The  order  of  these  two  bounding  0-D  domains  defines  a  direction  for 
which  is  used  to  set  the  direction  of  the  unit  tangent  vector  t  at  the  points  of 

the  arc.  The  corresponding  unit  normal  vectors  n  are  obtained  by  rotating  t 
counter-  clockwise  through  an  angle  n/ 2  (see  Figure  1.2.3). 


Figure  1.2.3 

The  particular  form  of  is  specified  by  its  signed  curvature  C^.  More  specif¬ 
ically,  if  Cl  =  0,  then  oj  is  a  straight  line  segment,  else  it  is  a  circular  arc  with 
radius  1/|C^  |.  In  the  latter  case,  if  C^>0  or  C^<0,  then  the  center  of  the  circle  is 

to  the  right  or  the  left  of  the  arc,  respectively,  when  looking  in  its  direction  (see 
Figure  1.2.3). 


Figure  1.2.4 


Hence  the  tangent  and  normal  vector  at  any  point  of  the  arc  are  given  by 


t 


cos  i!>  +  sin  d 
sin  i5  -  s^  cos  i3 


where 

a  x2 - xi  .  R  y2-y1 

cos  •&  =  — - —  ,  sm  £  =  — - 

L  =  J(x2-xl)2  +  (y2-yl)2  ,  p  =  yC 

S;  =  p  (  2  C  -  1)  ,  =  yj 1  -  p2(2C  -l)2 

Accordingly,  the  global  coordinates  (x,y)  and  the  local  coordinates  (£,ri  )  are 
related  as  follows 


_  - 

_  _ 

2^-1 

X(0 

i 

Xl+X2 

L 

cos  d  -sin  i3 

y®m 

“  2 

_yi  +  y2_ 

+  2 

sin  d  cos  d 

1  -  (2C-1) 

V  Vl-p2  +  yj  1  -  P2(2C  *  l)2 

Line  integration  of  a  function  f  along  the  1-D  domain  then  becomes 

Jfdy  =  lJ£<jc 

n1  0  ^ 

K 

Transformation  of  a  2-D  domain  is  defined  by  a  Coon's  patch;  that  is,  by  blend¬ 
ing  the  transformation  functions  (x(k},y(k))  ,  k=1,2,3,4  of  the  four  bounding  1-D 
domains.  As  Figure  1.2.5  shows,  the  0-D  and  1-D  domains  are  indexed  in  a 
counter  clockwise  direction. 


\ 


x=x(4,ti) 

y=y(4,ri) 


2  (2)  5 


Figure  1.2.5 

The  transformation  functions  can  be  expressed  in  the  following  form: 


-X1 

XV'(V 

\j 

-X4 

(2>(^) 

0 

x(4)(^4) 

0 

-X2 

x(3)(n3) 

0 

-X3 

y(^i)=[i-s  i  d  y()^2) 


0 

-yi 

y(1)0ii) 

0 

-y4 

y(2)(^2) 

0 

y(4)(^4) 

0 

(3).  . 

0 

-y2 

y  0i3) 

-y3 

Here  (xjO,yjo)i  j=i  ,2,3,4,  are  the  coordinates  of  the  four  0-D  domains  and 
(x(k)(rjk) ,y(k)(t|k)),  k=1,3,  (x(kHCk),yM(Ci<)).  k=2,4,  the  four  1*D  transformation 
functions.  Their  arguments  r|k,  or  C,k  reflect  the  appropriate  orientation  of  the 
corresponding  1-D  domains;  that  is, 


A 


I 

ft? 


1  k  *  l  k  *  ^ 

\  =lk^\+~2~~'>  >  ^k  =  lk^k+~2  — ) 


where  tk  =  +1,  if  the  orientation  of  the  k-th  1-D  domain  corresponds  to  the 
orientation  of  the  y-axis  or  x-axis,  respectively,  and  ik  =  -1  if  it  is  the  opposite. 
Since  the  determinant  of  the  Jacobian  of  this  transformation  appears  in  each 
volume  integral,  it  is  important  to  express  it  in  an  efficient  form  for  later 
evaluation.  The  partial  derivatives  can  be  written  as 


^  y(^,Ti) 


r2(4)  sin  i32 

x3(r()  -  XjCn) 

r4(4)  sin  134 

r2(5)  cos  i32 

y30i)  -  y*(Ti) 

r4(^)  cos  t34 

rj(Tl)  sin 

x4(^)  -  x2^) 

r3Cn)  sin  t33 

r/n)  cos 

y4(5)  -  y2(5) 

r3(ri)  cos  i33 

where 


rk(0  =  Lv 


Pk(2C*D 


1  -  pj  (2  ?  -  l)2 


'  Pk^)  =  Pk 


1  -  (2  C -1) 


1*P  l  +  J1-P2(2C-1)2 


^  m 

xk(C)  \  ^  COST^  -sin  ik  ( 2  C  -1) 


yt(C)  y[ 


%  COST^  Pk(Q 


0  0 

X.  +  X,  , 
k  k-1 


J  K  K-t 

Too 
yk  +  yk-i 


k  =  2,3,4  and 


o  o 
x,  +  x. 
1  1  4 


2  o  o 


1.3.  Finite  Element  Approximation 


The  mesh  on  our  domain  ft  consists  of  curvilinear  elements  which  are  first 
defined  on  the  basic  square  S  =  {  (^,ri)  ,  0  <  <  1}  and  then  mapped  onto  the 

closures  of  the  2-D  domains  which  make  up  ft.  A  nine-node  planar  Lagrangian 
quadratic  element  is  used  as  the  master  element  of  NFEARS. 


1(0,0) 

5  - 6 


Figure  1.3.1 

The  shape  functions  associated  with  its  nine  nodes  (see  Figure  1.3.1)  are  as 
follows: 


9i  =  ^MM)v(v-i) :  92  =  jC-^V^-i)  ;  (p3  =  ^u (n+i)v(v-i) 


95  =  0*M- 2)C''v’2)  ;  96  =  -^(i-v2)p  (h +i) 


1  1  2  1 
9-,  =  tM-  (M-  'I  )v(v+1);  (P  =-(1-p  )v(v+1);  <p -  -hi  (p  +1  )v(v+1 ) 


As  in  [26]  the  admissible  meshes  on  S  are  defined  recursively  by  the  two  rules  : 


(a)  The  mesh  consisting  of  the  four  quarter-squares 

Sk  =  {(^.h)  ;  i  ^  <  i+1 ,  j<2r|<j+1  },  k  =  t  +j+2i,  i,j=0,1 

of  S  is  admissible.  These  four  squares  S^Sg.Sg.S,*  are  called 
superelements. 

If  A  is  an  admissible  mesh  on  S,  then  the  mesh  is  admissible  that  is 
obtained  from  A  by  quartering  any  one  square  s  of  A  into  four  congruent 
squares  of  half  the  side-length  of  s. 


(b) 


S: 


Si 


v 


(5o  -  no) 


Figure  1.3.2 


Figure  1.3.2  shows  an  example  of  an  admissible  mesh  on  S.  Any  element  s  of 
an  admissible  mesh  A  on  S  is  uniquely  characterized  by  its  center  point  (^o.  ho) 
and  side-length  h  ,  and  the  transformation 


%  =  +  (0.5h)  p ,  v  =>  r|  =  rio  +  (0.5h)  v 
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maps  the  master  element  onto  s.  This  specifies  the  shape  functions  on  s. 
Instead  of  detailing  them,  we  give  here  the  more  general  formulas  for  quadratic 
interpolation  on  s.  As  above  let 

•  v  =  f  (1  "V 

and  suppose  that  Uj,  i=1,2 . 9  are  given  scalar  values  at  the  nine  nodes  of  s 

indexed  in  the  same  order  as  in  Figure  1.3.1.  Then  the  interpolated  value  u(^,r|) 
at  the  point  (S,,r|)  e  s  is  given  by 

■Ul-U2-U4+U5  2(U4-U5)  U7-Ug-U4+U5jr  y(v.1 } 

U(q,Tl)  =  [M-  (M-  ■"• ).  1.  M-  (p  +1  )]  2(U2'U5)  4u5  2(u8'u5^  1 

U3-U2'U6+U5  2(U6'U5)  U9-U8‘U6+U5  LV(V  +  1  } 

It  will  be  useful  to  record  also  the  corresponding  quadratic  interpolation  which  is 
induced  on  any  line  segment  [£o-h/2  ,  £o+h/2].  If  Uj.  i-1 ,2,3  are  the  values  at  the 
three  nodes  (indexed  in  order  of  increasing  ^-values),  then  the  value  u(£)  of  the 


one-dimensional  quadratic  interpolation  is  specified  by 


2 


u(C)  =[f(P-1).  1,  p(p+1)1  o 


) 


with  C,  =  q  or  r\. 


Once  the  shape  functions  are  defined  on  a  given  admissible  mesh  on  S,  we 

have  obtained,  at  the  same  time,  the  corresponding  space  of  conforming  finite 

element  functions  on  that  mesh.  Hence,  by  applying  the  transformations  from  S 

2 

onto  each  one  of  the  closed  2-D  domains  Ok  ,  we  construct  a  curvilinear  mesh 

on  the  full  domain  Q.  This  is  schematically  indicated  in  Figure  I..3  for  the  basic 
mesh  of  four  superelements. 
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Super-elements 


Figure  1.3.3 


Let  Vj  =  vj(x,y)  denote  the  global  shape  function  belonging  to  the  node  i  of  the 
resulting  mesh  on  Q.  Then  the  finite  element  space,  to  be  considered  here, 
consists  of  all  the  functions  of  the  form 


IN 

u(x,y)  =*  ^  UjV^x.y)  ,  x,ye  Q 


where  the  index  N  denotes  the  total  number  of  "free"  nodes;  that  is,  of  the  nodes 
not  contained  in  the  0-D  and  1-D  domains  carrying  Dirichlet  boundary 
condition. 


Dirichlet  conditions  can  be  prescribed  on  a  set  ©  consisting  of  0-D  and  1-D 
domains,  and  they  may  either  depend  linearly  on  a4  or  be  independent  of  that 
parameter;  that  is, 


.  .  /  b(x,y)  if(x,y)e  ©1 

u(x,y)  =  <  o4b(x,y)  if  (x,y)  e  ©2  @1  u  ©2  -  0 


Here  the  following  conditions  are  assumed  to  hold: 

(a)  The  function  b(x,y)  is  a  quadratic  polynomial  in  the  local  coordinate  £  on 
each  1-D  domain  in  A,  (see  Section  1.2  for  the  relevant  mapping) . 

(b)  The  intersection  0,  n  @2  may  contain  only  0-D  domains  with  zero 

prescribed  values,  u(x,y)  =  0 


(c)  If  a  1-D  domain  is  in  some  set  ©jt  (i=1,2),  then  the  two  adjacent  O-D 
domains  also  must  be  in  the  same  ©j. 

(d)  If  two  or  more  1-D  domains  with  a  common  adjacent  0-D  domain  are  in 
©j,  then  the  various  b(x,y)  functions  defined  on  these  1-D  domains  must 

have  an  identical  value  at  this  0-D  domain. 

The  functions  b(x,y),  as  quadratic  polynomials  in  their  local  coordinates,  are 
specified  by  their  values  at  the  0-D  domains,  and  at  the  mid-point  of  their  1-D 
domain  (see  Section  11.2).  The  necessary  quadratic  interpolation  for  all  other 
points  of  the  1-D  domain  is  performed  by  the  program. 


With  these  specifications  the  finite  element  approximation  of  our  given  problem, 
on  any  admissible  mesh,  is  a  system  of  N  nonlinear  equations: 


F 


(u)  =  J 


n 


v  +  iSL 

au  1  au 


-  |  G2  w  dxdy  -  J 


n 


3v: 


■v*,  9v, 
dO  i 


3x  3Uy  9y 


G1  Wj  dy  =  0 


dxdy 

,  i-1.2 . N 


where 


u  =  u(x,y)  is  the  desired  finite  element  approximation  that  satisfied  the 
boundary  conditions; 

v,  is  the  global  shape  function  corresponding  to  the  i-th  node  of  the 
(curvilinear)  mesh  on  Q; 


r  is  a  given  subset  of  the  union  of  0-D  domains  and  1-D  domains  on  which 
the  Neumann  boundary  conditions  specified  by  the  third  integral  are 
defined; 

Wj  is,  for  any  node  on  T,  the  induced  one-dimensional  quadratic  shape 
function; 


<J>  =  <D(0rx,y,U,Ux,Uy) ,  and  G.  =  Gj(aj,x,y)  ,  i=1 ,2 

are  the  given  functions  of  the  problem  as  discussed  in  Section  1.1 ; 


a*  ,  i=1 ,2,3  are  the  two-dimensional  parameter  vectors  and  04  the  scalar 
parameter  discussed  in  Section  1.1.  Recall  from  Section  1.1,  that  for 
the  computations  the  parameters  a-i,  02,  03  are  specified  as  linear 
functions  of  two  effective  parameters  X2. 

The  Jacobian  of  the  non-linear  mapping  F  =  (Fi,  F2 . Fn)  is  easily  obtained  by 

direct  differentiation.  For  ease  of  notation,  let 


0  1  3u  2  3u  0  1  ^vj  2  ^vi  .  _ 

u  =  u  ,  u  =  —  ,  u  =  —  ,  v  =v,v=  —  ,  v  =  —  ,i=1 ,2 N 

3x  3y  3x  3y 


Then  the  components  of  the  Jacobian  of  F  have  the  form 


r*  3o  k  1  ,  ,  . .  .  . 

L — r - r  V.  v  dxdy  ,  i,j=1 ,2 . N 


_  f  V  3  0  k  1 
aUj  =  J  ffi  3uk  3u'  Vi  Vj 


and,  with  the  notation  of  Section  1.1,  the  derivatives  with  respect  to  the  effective 
parameters  are 

3Fi  r  y[  Bi  32 0  +  ai  320  1  k 

1  n  K  UL  1  4  J 


i  f°'J 2  i 

-  pi,  J  -j  v,  dx  dy  -  PJ3J— L  w(  dy  (i=1 ,2,...,N;  j=1 ,2;  p>=p4,p2=0) 
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Here,  it  should  be  noted  that,  in  the  assembly  of  the  function  values  and  the 
Jacobian  components,  the  fixed  boundary  nodes  are  not  represented.  Thus  the 
derivatives  with  respect  to  o4  have  to  be  transfered  to  the  free  nodes. 


In  NFEARS  the  admissible  meshes  on  the  basic  square  S  =  {(£,r|),  0<^,r|<1}  are 
represented  by  a  combination  of  tree  structures.  This  representation  is  an 
extension  of  that  given  in  [25]  for  linear  elements.  More  specifically,  the 
quadratic  elements  of  the  mesh  in  the  2-D  domains  are  stored  as  extended 
quad-trees  and  their  one-dimensional  parts  on  the  free  1-D  domains  as  binary 
trees.  This  tree  structure  facilitates  the  required  mesh-refinements  and  de¬ 
refinements.  We  refer  to  the  cited  article  for  a  detailed  description  of  the  data 
structure  and  summarize  here  only  some  of  the  changes  that  were  needed  to 
accomodate  quadratic  elements. 

As  discussed  in  Section  1.3  above,  all  meshes  on  Q  are  first  defined  on  the 
basic  square  S.  Moreover,  by  definition  of  the  admissible  meshes,  the  initial 
mesh  on  each  2-D  domain  corresponds  to  a  basic  subdivision  of  S  into  four 
congruent  squares,  the  superelements.  Figure  1.4.1  shows  this  initial 
subdivision  with  the  four  parts  numbered  2  to  5.  The  mid-point  carries  the 
number  1  and  the  mid-points  of  the  sides  are  assigned  the  indices  6  to  9. 
Further  side  nodes  are  required  on  the  boundary  of  S  but  they  are  not 
numbered  in  the  Figure.  The  superelements  may  be  further  refined  but  can 
never  be  de-refined;  that  is,  each  2-D  domain  contains  at  least  nine  internal 
nodes. 

The  initial  sub-division  is  represented  by  a  labeled  tree,  where  the  root 
corresponds  to  the  mid-point  of  the  unit-square  and  the  8  descendant  nodes  to 
the  four  superelements  and  the  four  side-nodes,  respectively.  Labels  lx  and  ly 
with  values  -1 , 0  and  +1  are  assigned  to  the  arcs  of  the  tree  to  characterize  the 
geometrical  location  of  the  corresponding  node  in  relation  to  the  mid-point.  Note 
that  the  terminal  nodes  2  to  5  represent  mid-points  of  open  squares,  nodes  6  to 
9  are  mid-points  of  open  intervals,  and  the  root  corresponds  to  a  single  point. 


Subdivided  2-D  subdomain: 


Corresponding  tree: 


The  thin  lines  in  Figure  1.4.1  indicate  the  effect  of  a  refinement  of  elements  2  and 
4.  Eight  descendant  nodes  are  appended  to  the  particular  tree  nodes  in 
correspondance  with  the  four  new  elements  and  four  new  side-nodes.  As  in  the 
case  of  the  arcs  from  node  1 ,  the  eight  arcs  to  the  new  nodes  carry  labels  lx,  ly 
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with  values  -1.0.+1.  Two  further  side-nodes  are  introduced  on  both  sides  of 
node  8.  These  are  implemented  by  adding  two  son-nodes  to  node  8  and 
labeling  the  arcs  |x=0,ly=-1  (below  8)  and  |x=o,ly=+1  (above  8),  respectively.  In 
addition,  so-called  irregular  nodes  appear  on  both  sides  of  nodes  6  and  7. 
Generally,  a  node  of  any  mesh  is  irregular  if  it  lies  on  the  side  of  some  element 
but  is  not  the  mid-point  of  that  side.  These  irregular  nodes  are  not  represented 
in  the  tree.  Solution  values  at  the  irregular  points  are  linearly  dependent  on  the 
solution  values  at  the  corner-points  and  at  the  mid-point  of  the  side.  Side  points 
are  also  introduced  on  the  bordering  1-D  domains.  Whether  these  are  regular 
nodes  or  irregular  nodes  will  depend  on  the  mesh  on  the  neighbouring  2-D 
domain,  if  there  is  one.  If  there  is  none;  that  is,  if  the  1-D  domain  is  an  outside 
boundary,  then  those  nodes  are  defined  to  be  regular  nodes. 


The  resulting  labeled  tree  has  two  important  properties: 


1.  All  non-terminal  nodes  of  the  tree  represent  single  points  while  the 
terminal  nodes  correspond  to  the  mid-points  of  open  squares  or  open  horizontal 
or  vertical  intervals  in  the  local  coordinate  system.  Open  squares  are  mapped 
onto  the  interior  of  elements  of  Q  and  the  lx.ly  labels  on  their  incoming  arcs  carry 
non-zero  values.  Open  intervals  are  those  which  have  a  label  lx=0  (vertical 
sides)  or  lv=0  (horizontal  sides). 

2.  There  exists  a  unique  path  from  the  root  to  any  given  node  on  the  tree 
which  consists  of  arcs  with  labels  Uj.ly,  ,  i=1,2,...,m,  that  define  the  local 

coordinates  of  the  corresponding  point  as 


m  m 

5-7<1+Z  '‘i2")  and  l-xO+X  lv,2") 


For  terminal  nodes  the  length  m  of  the  path  specifies  the  length  h  of  the  side  of 


the  corresponding  open  square  or  interval  as  h  =  2 


The  1-D  domains  are  mapped  onto  open  unit  intervals  (see  Section  1.3).  Each 
free  1-D  domain  is  represented  by  a  labeled  binary  tree,  while  fixed  1-D 
domains;  that  is,  those  that  carry  Dirichlet  conditions,  are  only  documented  in 
summary  records.  Initially,  a  binary  tree  has  one  root-node  and  two  son-nodes 
with  arcs  carrying  the  labels  -1  and  +1,  respectively.  The  root  node  corresponds 


w  -r.  •f.  a  «"» A.A. 


to  the  mid-point  of  the  1-D  domain  and  is  also  a  cornerpoint  of  a  superelement 
in  the  adjacent  2-D  domains;  the  son-nodes  are  then  representing  the  sides  of 
this  superelement.  This  is  shown  in  Figure  1.4.2. 


1-0  subdomain  separating  two  2-D  subdomains: 


O'D  subdomain 


V 


i 

v. 


N" 


I.- 


For  the  control  of  the  mesh  modification  process,  as  well  as  for  the  initial 
definition  of  the  meshes  on  £3,  we  follow  [3],  and  introduce  a  density  function 
and  mesh  intensity. 

2 

The  2-D  domain  Qk  is  mapped  onto  a  unit  square  which  contains  four 
superelements  S1*,  j=1 ,2,3,4.  As  in  Section  1.3  let  Vj(£,r|),  i=1,...,9,  be  the  nine 

k 

global  shape  functions  of  the  superelement  Sr  For  every  superelement  exactly 

one  corner  point  is  the  image  of  a  0-D  domain.  Suppose  that  this  point 
corresponds  to  the  corner  with  index  I  of  the  master  element.  Moreover,  let 

k 

ro(q,n)  denote  the  distance  between  it  and  any  point  (^,n)  e  S,.  Then  the  un- 

k  k 

normalized  density  function  dj  on  Sj  is  defined  by 

9 

log(d|'(4.n))  =  X  p,  v<5’1)  *  2  Pa  v|(vn)  tog(r0(^,n»  .  i.n  e  s* 


i 

£ 

£ 

I 


where  po,...,p9  are  given  constants.  This  definition  pre-supposes  that  the 
solution  can  be  split  into  a  smooth  part  and  the  effect  of  a  singularity  near  the  0- 
D  domain  associated  with  S*  .  More  specifically,  the  sum  over  the  shape- 

functions  reflects  mesh  for  the  smooth  part  and  the  logarithmic  term  mesh  for 
the  singular  part.  Line  singularities  are  not  represented  separately;  that  is,  in 
our  setting,  they  must  be  incorporated  into  the  smooth  part. 


B 

;8 

l 


ft 


f 

fc» 


When  a  node  belongs  to  two  different  superelements,  then  the  coefficients  pj 
corresponding  to  it  in  these  elements  are  assumed  to  be  identical.  With  this  a 
continuous,  positive  density  function  dk(^,ri)  is  defined  on  the  entire  unit  square. 
Moreover,  dk(£,,n)  is  jniquely  characterized  by  29  constants,  namely,  25 

coefficients  pi . pg  (nine  on  the  open  2-D  domain,  three  each  on  its  four  open 

1-D  domains,  and  one  each  at  its  four  0-D  domains  ),  and  4  singularity 
exponents  po  (one  each  at  the  four  0-D  domains). 

The  density  functions  dk  can  be  extended  to  a  continuous,  positive  density 
function  d  on  the  entire  domain  Q  by  requiring  that  the  coefficients  at  any  nodes 
m  the  intersection  of  two  closed  2-D  domains  are  identical.  This  function  d  is 
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completely  characterized  by  nine  coefficients  for  each  open  2-D  domain,  three 
for  each  open  1-D  domain,  and  two  for  each  0-D  domain.  Then 


Dn  (M  = 


d  (^,r|) 


shall  be  the  normalized  density  function  on  the  domain  Cl. 

Let  Aj,  j  =  1 ,2,...,  be  a  sequence  of  meshes  with  Ej  elements  such  that  Ej  -* «  as 
j  oo.  Moreover,  assume  that  for  any  open  subset  O'  of  Q  which  contains  EjU. 
(open)  elements,  we  have  the  inequality 


where 


°S (j.^)  J DnM)  drj  <  -^p-  <  a2(j,ft')  J  D^^.ri)  d^  dri 


n 


n- 


lim  sup  a2(j,Q)  and  lim  inf  a  (j,Q) 

i  j  ->  oo 


are  finite  and  independent  of  Cl'.  Then,  evidently,  for  large  j  the  quantity 

3(ff)  -  A-  J  Dn(5,n)  d5  dn 

iQ'  n' 

is  nearly  independent  of  Q'.  This  suggests  that  we  may  characterize  a  mesh  by 
its  normalized  density  function  and  a  number  3  ,  called  its  intensity,  such  that 

En-  =4 1  Do«-n)<*s«m 

a 


is  the  number  of  elements  in  any  open  subset  ft'  of  Cl.  This  relation  allows  us  to 
construct  an  admissible  mesh  A  on  Cl  once  a  normalized  density  function 

DQ(£,n)  and  an  intensity  3  are  known.  Clearly,  it  suffices  to  give  this  algorithm 


2  2 

only  for  some  closed  2-D  domain  Qk  of  Q  .  As  before,  let  £2k  be  mapped  onto 

the  unit  square  Sk  and  let  Ak  be  the  mesh  consisting  of  the  four  superelements 
Sk. 

For  all  elements  co  of  Ak  dfi 

Compute  for  the  set  £1  =  {co}; 
tf  Eft-  >  1  then 

construct  a  new  mesh  Ak  by  subdividing  co  into  four  new 
elements, 

else 

terminate  the  process  with  the  current  mesh  Ak. 

The  process  ends  with  a  mesh  Ak  for  which  En.  <  1  for  all  singleton  sets  £T  =  {co} 
of  Ak.  Moreover,  we  also  know  that  En,  >  1  for  the  father  co'  of  co;  that  is,  for  the 
element  of  an  earlier  mesh  from  which  co  was  constructed  by  the  quartering 
rule. 

Some  examples  of  this  process  are  given  in  Figure  11.3.4  .  It  may  be  noted  that 
there  the  intensity  was  set  equal  to  1/20,  and  hence  we  should  expect  a  mesh 
with  20  elements.  But  the  actual  number  of  elements  is  closer  to  50  due  to  our 
use  of  the  quartering  rule;  that  is,  the  requirement  that  only  admissible  meshes 
are  allowed. 


In  Sections  1.3  we  specified  our  problem  in  the  form  of  a  system  of  nonlinear 
equations 

Ffu.^.J^)  =0  ,  F:  RN  x  R2  ->  RN, 

involving  the  vector  u  of  the  unknown  values  at  the  N  free  nodes  of  the  finite 
element  mesh  on  Q  and  the  two  effective  parameters  X1 ,  X2.  This  system 

approximates  the  nonlinear  operator  equation  introduced  in  Section  1.1,  and,  as 
indicated  there,  the  set  of  all  solutions  (u,\i,X2)  e  RN  x  R2  forms  a  differentiable 
manifold  in  RN  x  R2.  We  summarize  here  briefly  some  of  the  necessary 
differential  geometric  properties  of  this  solution  manifold  and  refer  for  further 
details,  for  instance,  to  [19]  and  the  literature  cited  there. 

For  ease  of  notation,  we  shall  use  the  abbreviation  x  =  (u,  XlP  X2)  e  RN+2  and 
write  DF(x)  for  the  N  x  (N+2)  Jacobian  of  F.  A  point  x  is  called  regular  if  the 
Jacobian  has  full  rank  N.  Then  the  set  of  all  regular  solutions 

M  =  { *  e  RN+2;  F(x)  «  0,  rgeDF(x)  =  N  } 

is  a  two-dimensional  differentiable  manifold  in  RN+2  of  the  same  differentiability 
class  as  F.  By  restricting  ourselves  to  this  regular  solution  manifold  we  asume 
that  a  suitable  unfolding  of  the  original  problem  has  been  chosen.  In  practical 
applications,  this  is,  in  general,  a  natural  assumption. 

As  indicated  already  in  Section  1.1,  we  wish  to  analyse  the  shape  and 
characteristic  features  of  this  solution  manifold.  However,  for  this  we  have  to 
recall  that  the  mapping  F  represents  only  a  discrete  approximation  of  the 
original  (infinite)  dimensional  problem  of  Section  1.1.  As  noted  there,  this 
original  problem  --  after  restriction  to  the  two  effective  parameters  Xi,  X2  -  may 
also  be  expected  to  possess  a  two-dimensional  solution  manifold  in  the  product 
of  the  trial  space  X  and  the  parameter  space  R2.  It  is  actually  this  original 
manifold  which  we  wish  to  analyse;  but,  of  course,  that  manifold  is  not  directly 
accessible  for  a  computational  study.  This  raises  the  question  of  the  influence  of 
the  finite  element  discretization  upon  the  shape  and  features  of  the  original 
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manifold  and  of  estimates  of  the  size  of  the  discretization  error.  This  will  be 
discussed  in  Section  1.8  below.  Here  we  concentrate  on  the  above  manifold  M 
defined  by  the  specific  finite  element  discretization  F. 

At  any  point  *  e  M  the  tangent  space  T*M  of  M  may  be  identified  with  the  null- 
space  of  the  Jacobian;  that  is, 

TXM  =  Ker  DF(x)  =  {  w  e  RN+2,  DF(x)w  =  0  }. 

Accordingly,  the  orthogonal  complement 

Nx  M=  (Tx_M)^  =  rge  DF(x)T 

is  the  normal  space  at  that  point. 

For  the  analysis  of  any  differentiable  manifold  we  require  suitable  local 
coordinate  systems  on  M.  For  the  computations  we  use  local  coordinate 
systems  induced  by  given  2-dimensional  subspaces  T  of  RN+2.  More 
specifically,  any  such  space  T  induces  a  local  coordinate  system  of  M  at  x  e  M 
provided  that 

To  TxM  ={0}. 

If  this  condition  holds  then  there  exist  open  neighborhoods  V1  and  V2  of  the 

origins  of  T  and  RN+2  ,  respectively,  as  well  as  a  unique  differentiable  function  w: 
Vi  ->  T1  ,  w(0)  =0,  such  that 

M n V2  =  {  ye  RN+2  ;  y  =  x  + 1  +  w(t),  t  e  V1 }. 

A  point  x  is  a  non-singular  point  with  respect  to  the  given  coordinate  space  T  if 
T  induces  a  local  coordinate  system  of  M  at  x  ,  otherwise  it  is  called  a  singular 
point  or  foldpoint. 

Let  A  be  the  two-dimensional  subspace  of  RN+2  spanned  by  the  effective 
parameter  directions.  Then  interest  centers  especially  on  determining  the 
foldpoints  with  respect  to  the  coordinate  space  T  =  A  .  If  our  original  problem 

27 


,»r,*»vrPr^r4*l 


WWW 


concerns  a  study  of  equilibrium  configurations  of  a  mechanical  structure  then 
these  foldpoints  with  respect  to  A  represent  approximations  to  the  eqilibria 
where  a  change  of  the  stability  behavior  of  the  structure  may  be  expected  occur. 

The  basic  procedures  for  the  computational  analysis  of  our  solution  manifold  M 
are  the  continuation  methods.  These  methods  require  a  restriction  to  some  path 
on  M  and  then  produce  a  sequence  of  points  along  that  path.  One  such 
continuation  procedure  is  incorporated  in  NFEARS.  The  allowable  paths  for  this 
procedure  are  defined  by  specifying  the  effective  parameters  Xi,  X.2  as  linear 
functions  of  a  scalar  parameter  X: 

\  =  8j  X  ,  i=1 ,2  ,  51  +  S2  =1. 

Some  of  the  details  of  the  continuation  process  and  of  its  capabilities  are 
discussed  in  Section  1.6. 

Obviously,  it  is  not  easy  to  develop  a  good  picture  of  a  two-dimensional 
manifold  solely  from  information  along  such  paths.  This  led  recently  to  the 
development  of  methods  for  the  computation  of  simplicial  approximations  of 
two-dimensional  open  subsets  of  M.  NFEARS  incorporates  a  form  of  the  method 
introduced  in  [21], [22]  which  will  be  discussed  in  Section  1.7  below. 


1.6.  The  Continuation  Process 


The  continuation  method  incorporated  into  NFEARS  follows  the  design  of  the 
PITCON  system  (see  [23], [24])  .  We  summarize  here  only  the  principal  aspects 
of  that  method  and  refer  to  the  original  articles  for  further  detail. 

As  discussed  in  Section  1.3,  our  finite-element  approximation  has  the  form 


Fi(u,Xi,X2)  =  0  ,  i=1 .2....N,  ueRN 

Let  M  again  be  the  corresponding  regular  solution  manifold  and,  for  ease  of 
notation,  write  x  =  (u,  Xi,  X2).  For  the  continuation  process  we  have  to  restrict 

ourselves  to  a  path  on  M  through  a  specified  point  xP  =  (u0,X.°,X°)  e  M.  This  is 

equivalent  with  an  augmentation  of  the  system  by  some  scalar  equation.  As 
mentioned  in  Section  1.5,  NFEARS  uses  only  paths  on  M  which  are  specified  by 
a  linear  combination  of  the  effective  parameters  \-\,  ?i2.  In  other  words,  we 
consider  augmented  systems  of  the  form 

[  1 


f*(uA,.x.2)  = 


Fn(u,A.i,X>2) 
S1(Xl-X°)  +52(X2-X°) 


where  Si,  S2  are  the  given  constants  introduced  at  the  end  of  Section  1.5.  Let 

n  =  {  x  e  RN+2,  F*(x)  =  0,  rge  DF*(x)  =  RN+1  } 

be  the  regular  solution  manifold  of  the  augmented  equations.  The  connected 
component  no  of  n  containing  the  given  initial  point  x°  is  the  path  on  M  which  is 
to  be  computed. 

The  continuation  process  begins  from  the  given  starting  point  and  produces 
a  sequence  of  points  xk  ,  k=0,1 . which  approximate  points  of  M.  For  any  k  >  0 


the  step  from  xk  to  xk+1  corresponds  to  an  implementation  of  the  local 
coordinate  representation  discussed  in  Section  i.4.  More  specifically,  if 
T=span{t},  with  t  e  RN+2  ,  t  *  0,  is  a  local  coordinate  space  of  M  at  xk,  then,  for 
any  fixed  y  e  RN+2,  the  Jacobian  of  the  augmented  equations 

F*U) 

tV_*) 

is  non-singular  for  all  &  in  some  neighborhood  of  &k.  Thus  if  y  approximates  a 
point  of  n  in  that  neighborhood,  then  it  follows  readily  that  the  above  system  has 
a  unique  solution  xk+1  e  M  which  can  be  computed  by  means  of  a  locally 
convergent  iterative  process,  started,  say,  at  y  . 

This  gives  an  outline  of  the  process  and  it  remains  only  to  specify  the  particular 
choices  used  in  PITCON  and  NFEARS.  At  *k  the  normalized  tangent  vector  &k 
of  n  is  computed;  that  is, 

DF*(*V=  0  ,  lrt=1, 

where  for  k=0  the  direction  of  jyk  is  user-given  and  for  k>0  it  is  defined  by 
(  ^  )T  (*k  -  2ik'1)  >  0. 

As  a  check  against  leaving  the  connected  component  n0  the  condition 
det(DF‘(*k)T,  yyk)  >0 

is  also  monitored.  The  predicted  point  y  is  now  computed  by  linear 
extrapolation  along  the  tangent  direction;  that  is, 

y  =  £k  +  x  wK 

The  choice  of  the  step  length  x  is  the  same  as  in  PITCON  and  we  refer  to  [23]  for 
the  detailed  discussion. 


For  the  definition  of  the  local  coordinate  system  T  =  span{t}  at  xk  one  of  the 
natural  basis  vectors  of  RN+2  is  used;  that  is  we  set  t  =  §}  where  the  index  i  is 
chosen  such  that  |(wk)T£'|  is  maximal.  This  defines  the  above  given  augmented 
system  which  is  then  solved  by  means  of  a  chord  Newton  process  started  from 
the  predicted  point  y.  For  details  of  the  acceptance  and  rejection  tests  of  this 
iterative  process  and  for  the  recovery  procedures  after  failure  we  again  refer  to 
the  cited  article  [23]. 

The  process  allows  for  the  computation  of  target  points;  that  is  of  points  x  e  M 
where  the  component  xTei  with  a  prescribed  index  j  has  a  specified  value 
In  order  to  detect  such  a  point  the  condition 

sign(  (x  k-^Tei  -  )  *  sign(  (x  k)Tei  -  ) 

is  monitored.  If  it  holds,  the  interpolated  point 

y'rO-uJxk-Uuxk  ,  (ei)T y*  =  g* 

is  evaluated  and  then  the  corrector  process  is  applied  with  the  augmenting 
equation  (ei)Tx  -  =  0  and  with  y*  as  starting  point.  Once  again  we  refer  to  [23] 

for  further  details. 

Finally,  the  continuation  process  incorporates  a  procedure  for  computing  limit 
points  of  n  with  respect  to  a  specified  natural  basis  direction  ei;  that  is  of  points  x 
on  n  where  the  j-th  component  of  the  tangent  vector  is  zero.  The  computation  of 
these  points  is  similar  to  that  of  the  target  points.  For  their  detection  the 
condition 

sign(  (wk-i  )Tei)  *  sign(  (wk)Tei) 

is  monitored.  If  it  holds,  then  the  limit  point  process  of  [23]  in  the  form  discussed 
in  [15]  is  applied.  We  refer  to  these  papers  for  the  details. 


1.7.  Simolicial  Approximations  of  the  Solution  Manifold 


As  observed  already  earlier,  it  is  difficult  to  gain  a  picture  of  the  two-dimensional 
solution  manifold  M  of  our  nonlinear  equations 

F(x)  =  0  ,  x  e  RN+2 

solely  from  the  information  along  the  paths  used  by  the  continuation  process. 
For  this  purpose  a  form  of  the  method  described  in  [21],  [22]  for  computing 
simplicial  approximations  of  open  regions  of  M  was  incorporated  in  NFEARS. 
Once  again  we  refer  to  the  original  articles  for  details  and  concentrate  mainiy  on 
the  features  that  are  different. 

The  approximation  process  begins  with  the  choice  of  a  reference  mesh  A  of  R2. 
In  the  original  method  a  equilateral  triangular  mesh  was  used  while  NFEARS 
employs  a  planar  Kuhn-triangulation  A  as  shown  in  Figure  1.7.1. 


Figure  1.7.1 


The  aim  of  the  triangulation  process  is  to  transfer  the  nodes  of  a  part  of  A, 
together  with  their  connnectivity  information,  from  R2  onto  M.  This  transfer 
process  works  with  groups  of  nodes.  More  specifically,  in  NFEARS,  we  use 
"patches"  of  eight  triangles  which  form  rectangles  of  A  and  contain  nine  nodes, 


namely,  one  center  node  and  eight  boundary  nodes.  Figure  11.6.1  shows  a 
portion  of  such  a  triangulation  with  one  patch  of  triangles  indicated  as  a  hatched 
rectangle. 

The  continuation  process  is  used  to  generate  a  first  "center"  point  x  e  M  for  the 
"region"  on  M  that  is  to  be  triangulated.  Now  a  patch  of  the  reference  mesh  A  is 
mapped  linearly  onto  the  affine  tangent  space  x+T^M  in  such  a  way  that  the 
center  point  of  the  patch  corresponds  to  x.  For  this  mapping  we  use  an 
appropriate  basis  w-i,  W2  of  TXM  together  with  two  user-defined  scale  factors  xi, 
t2  for  each  one  of  these  directions.  In  other  words,  the  images  of  the  nine  nodes 
of  the  patch  on  x+TxM  are 

x  +  ki  Tiwi  +  k2X2W2  ,  ki=  -1 ,0,1  ,  k2=-1 ,0,1 
These  points  are  now  projected  onto  M.  In  the  original  method  in  [21]  ,  [22]  a 
chord  Gauss-Newton  method  is  used  to  produce  a  projection  from  x+T^M  onto 
M  which  is  orthogonal  to  the  tangent  space.  In  order  to  retain  the  same  data 
structure  as  in  the  continuation  process,  we  use  in  NFEARS  the  chord-Newton 
process  applied  to  the  augmented  equations 

F(Z) 

(e')T(z  -  y)  =0 

.(e])T(z-  y)_ 

Here  i  and  j  are  the  indices  of  the  largest  components  in  modulus  of  w1  and  w2, 
respectively,  and  y  e  x+T^M  is  the  point  which  is  to  be  "projected"  onto  M. 

In  order  to  repeat  these  steps  with  a  neighboring  patch,  a  new  center  point  on  M 
is  constructed  by  projecting  an  appropriate  point  along  one  of  the  basis  lines  of 
i+T^M  onto  M.  More  specifically,  if,  in  A,  the  new  patch  is  located  in  the 
"upward",  "left",  "right", or  "downward"  direction  from  the  original  one,  then  the 
new  center  point  is  the  projection  onto  M  of  the  points  &+2x2W2,  x-2x-|Wi, 
X+2t-|Wi,  *-2t2W2,  respectively.  Now  a  suitable  basis  of  the  tangent  space  of  M 
at  the  new  center  point  has  to  be  constructed  which  then  allows  us  to  map  the 
new  patch  onto  M  (using  the  same  scale  factors  xi,X2  as  before).  But,  of 
course,  nodes  of  A  that  have  already  been  mapped  onto  M  will  not  be  used 
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Clearly,  in  order  to  ensure  that  the  images  of  the  rectangles  of  A  on  M  form  a 
curvi-linear  mesh  on  the  manifold,  we  need  to  use  bases  on  the  tangent  spaces 
TXM  which  change  smoothly  from  point  to  point.  In  the  terminology  of  differential 
geometry  this  means  that  we  require  an  algorithm  which  generates  a  moving 
frame  on  the  portion  of  M  under  consideration.  Standard  methods  for  computing 
tangent  bases  do  not  produce  continuously  varying  tangent  bases  (see  e.g. 
[11]).  A  slightly  modified  form  of  the  moving  frame  algorithm  introduced  in  [21]  is 
used  in  NFEARS  to  construct  the  basis  vectors  of  the  tangent  spaces  at  the 
center  points  of  the  patches. 


Suppose  that  x  e  M  is  a  point  where  a  basis  of  the  tangent  space  TXM  is  to  be 
computed.  If  we  choose  two  distinct  natural  basis  vector  e'  and  ei  of  RN+2, 
neither  of  which  is  orthogonal  to  T*M,  then  the  solutions  q1t  q2  of  the 

augmented  system 

r DFooi 


(e) 


N+k 

e 


k=1 ,2 


certainly  form  a  basis  of  T^M.  Accordingly,  the  vectors  t1  ,t2 


<?i 

II  qn ll2 


((q/q,  )qj, 


X  =  II  q2  -  ((q2>  q,)  ^\\2 


represent  an  orthogonal  basis  of  T^M.  As  in  the  continuation  process,  the 
indices  of  the  largest  components  in  modulus  of  t1  and  t2  can  be  used  to 

determine  the  indices  i  and  j  of  the  basis  vectors  in  the  augmented  system. 


Once  an  orthogonal  basis  t1  ,  t2  of  TXM  has  been  obtained  the  quantities 
0k  =  cos(dk)  =  ^((eV  t/  +  ((ek)Tt2)2  ,  k=  1,2 . N+2 
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are  the  direction  cosines  of  the  principal  angles  between  T^M  and  the  natural 
basis  vector  of  RN+2.  Let  ii,  12  be  the  indices  corresponding  to  the  largest  of 
these  0k,  with  ties  broken  by  lexicographic  ordering.  Then  we  form  the  matrix 


,  > 


ai1  ai 2 


a21  a22 


,  a,  =  (e‘ )  tk  ,  j.k  =  1,2 


and  with  it  the  quantities 


1  1 


2  1 


^  d  ^11  +  a22^  ’  5C  d  (ai2  ‘  a21 Z  -  Ay  (an  +  a22^  +  ^ai2  "  32V 


It  was  shown  in  [22],  that  the  vectors  w,,w2 


w,  =±(Z1  t1  -Z2t2)  ,  w2  =  ±  (x2  t1  +  x1  t2) 


with  the  signs  chosen  such  that 


(e'k)Ttk  >  0  ,  k=1 ,2 


provide  the  desired  "moving  frame"  basis  of  TXM.  Evidently,  only  the  original 
basis  vectors  qi  ,  q2  together  with  the  various  transformation  coefficients  have 
to  be  retained  and  the  two  sets  of  w-vectors  need  never  be  computed  explicitly. 


The  sequencing  of  the  "transfer"  of  the  patches  of  A  onto  M  is  controlled  by  the 
user  and  will  be  described  in  Section  11.6. 


Let  A  be  a  given  mesh  on  Q  and  u  the  computed  finite  element  solution  on  A. 
The  calculation  of  the  a  posteriori  error  estimate  e(A)  of  u  follows  the  approach 
for  linear  problems  discussed,  for  instance,  in  the  surveys  [2],  [7],  [8],  [10]  and  in 
[3],  [9].  In  particular,  e(A)  is  obtained  as  a  sum  of  squares  of  error  indicators  p(co) 
of  the  elements  co  of  A;  that  is, 

e(A)2  =  X  P(°*2 

These  error  indicators  are  calculated  in  a  two-step  process.  First,  it  is  assumed 
that  the  solution  is  sufficiently  smooth  and  hence  that  the  indicators  of 
neighboring  elements  do  not  differ  too  much.  Then  this  assumption  is  checked 
and,  if  the  indicators  of  neighboring  elements  (of  a  common  father  node  in  the 
quad-tree  of  Section  1.4)  are  different,  then  their  computed  indicators  are 
adjusted  accordingly.  This  adjustment  turns  out  to  be  desirable  since  the 
unmodified  error  indicators  under-estimate  the  error  near  singularities. 


In  addition  to  the  quadratic  shape  functions  of  Section  1.3  we  introduce  on  the 
master  element  of  Figure  1.3.1  the  error  shape  functions 

y/p.v)  =  v  (v2-  1)  (p2- 1)  ,  v|/2(p,v)  =  p  (v2  -  1)  (p2-  1) 

On  any  element  co,  this  induces  the  global  error  shape  functions,  W|(x,y)  = 
¥i(p(x,y),v(x,y)),  i=1 ,2,  and  with 


dX  ^ 

Rk  =  “A- - ,  k.1,2 

dx  dy 


Qk  =  J  (  X  u'  wk  '  G2(a2,x,y)  wk)  dx  dy  ,  k=1,2 

co  '.j=0  du  du  ‘ 


the  first  approximation  of  the  error  indicator  on  co  is  defined  by 


r  j 


JV 
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P'(0)  )  =  720 


Q?  aj_ 

r;  +  r= 


In  addition,  NFEARS  calculates  a  linear  energy  term  for  a>  in  the  form 


En  (co)  =  J  ^  ^  O  .  u  1  u  1  (jx  dy 

a,  i.j=°  3u  9u  1 


Now,  an  adjusted  error  indicator  p  =  p (co)  is  calculated  from  p'  for  those 
elements  co  whose  "brothers"  are  also  elements  of  the  current  mesh  A.  In  other 
words,  if  co'  is  the  father  of  co  then  all  four  "sons"  of  co'  must  also  be  elements  of 
A.  For  instance,  in  Figure  1.8.1,  only  the  elements  co,,  i=1 ,2,3  and  4,  are 
candidates  for  adjustment,  while  the  elements  marked  by  co0  retain  their  already 
calculated  error  indicators  p';  that  is,  we  set  p(coo)  =  p’(o>o)- 


Figure  1.8.1 

Any  such  set  of  four  elements  coi . 004  of  A  will  be  called  an  adjustable  cluster. 

The  adjustment  of  the  error  indicators  of  such  a  cluster  proceeds  as  follows:  We 
order  the  four  elements  in  decreasing  order  of  their  calculated  indicators  pi'  = 
p'(coj)  ,  say, 
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Pi'  ~  P2'  ~  P3'  -  P4 

and  then  define  the  ratio  between  the  two  largest  error  indicators  and  those  of 
the  corresponding,  diagonally  opposite  elements  co^  and  co^: 


Pi’ 

P’(f>i) 


and 


p£ 

'2  P'(O)J) 


Now  the  adjusted  error  indicators  p  are  calculated  by  the  algorithm: 


if  P1'>8p2’ 

th£Q  p,  =  g-,(r1 ),  pk  =  pk'  ,k=2,3,4 
else  if  p2'  >  100  p3' 

then  Pi  =  g2(ri)*  P2  =  92(r2)'  P3  =  P3-  P4  =  P4' 

else  pk  =  pk'  ,k=l  ,2,3,4 

Here  the  two  function  gi(r)  and  g2(r)  are  defined  by 


n  _  /  3.481 3  e0  003762r  if  r  >  234 
9  i(r;  ”  \0.031 605  r  +  1  if  r  <  234 


{.  _ .  0.001  r  ___ 

1.7132  e  if  r  >  700 

0.0034999  r  +  1  if  r  <  700 


Let  h=h(w)  be  the  side-length  of  the  image  of  the  element  to  in  the  unit-square 
S.  Then  the  computed  error  indicator  p(co)  of  to  can  be  shown  to  have  the 
following  relation  to  the  H3(co)-norm  of  the  exact  solution: 

p2(fo)  =  C||u|f3  h4(<o) 

H  (<B) 


This  suggests  that  we  introduce  the  quantity 


x(co)  = 


J2I&L 

h(co)3 


for  which 
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in  the  local  coordinates  {^,r\)  on  the  unit  square  S. 


Let  Sk  be  the  unit  square  containing  the  element  co  and  suppose  that  a 
normalized  density  function  Dn(^,r])  and  intensity  3  are  given.  Then  the  size  h  of 

an  element  co  containing  the  point  is  asymptotically  equal  to 


h(co)  = 


.1/2 


D>.i> 


which  leads  to  the  asymptotic  error-indicator  formula 

t(co)2 


p"(co)2=  32J^f  d^dr!  .  e(A)2  =  X  P"^)2 

<o  Dn  “ e  A 

For  a  minimal  error,  this  formula  suggests  that  we  should  choose  the  density 
function  DQ(^,r|)  so  as  to  minimize  the  expression 

N 


I J 

k-i  i 


TftP) 

D? 


d^dti 


n 


This  minimization  is  performed  on  the  superelement  Sk  to  which  co  belongs.  The 

k  k 

definition  of  the  unnormalized  density  function  dj  on  a  superelement  Sj  is 

given  in  Section  1.4.  The  desired  coefficients  p0  and  pj  ,  i=1 . 9  are  calculated 

by  a  least  square  fit.  More  specifically,  for  any  element  co  in  Sj  with  side-length 

h(co)  we  introduce  the  least  squares  functional 
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which  is  minimized  to  calculate  p0  and  pi,...,p9. 


The  resulting  (raw)  density  functions  d*  over  the  superlements  of  Q  are  then 

adjusted  to  produce  a  continuous  function  over  the  full  domain.  If  a  node  occurs 
as  point  i,  (1  <i<9),  in  one  superelement  Sk  and  as  point  j,  (1  <j<9),  in  another 
superelement  Sk‘,  and  if  pj(co)  and  pj(oo')  are  the  corresponding  coefficients  of 


the  density  in  co  and  co',  respectively,  then  the  values  of  these  two  coefficients 
are  replaced  by  their  maximum;  that  is,  we  set 

Pi(Sk)  :=  Pj(Sk')  :=  max  {  pmin  ,  Pi(Sk),  pj(Sk') } 

with  a  suitable  minimal  value  pmin  •  Similarly  we  adjust  the  coefficients  p0  = 
p0(Sk),  but  now  we  use  minima.  In  other  words,  if  a  particular  0-D  domain  is  a 
corner  point  of  the  superelements  Sk  and  Sk'  then  we  set 


p0(Sk)  :=  min  {  pmax,  p0(Sk),  p0(Sk') } 


These  adjusted  p0 -coefficients  are  associated  with  the  four  0-D  domains  of  S 
and  with  this  the  calculation  of  the  unnormalized  density  function  dn(^,ri)  is 

complete.  The  un-normalized  density  function  is  normalized  again  as  described 
in  Section  1.4: 

dn  (ih) 

r'v  /fc  _ \  ** 


dn  (S,h)  drj 


With  this  adjustment,  each  closed  2-D  domain  has,  as  required,  exactly  29 
associated  coefficients  (see  Section  1.4).  Accordingly,  as  discussed  above,  we 
can  calculate  on  Q  the  error  indicator 

n. 
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where  the  average  intensity  3  should  be  chosen  such  that  en  <  emax  with  an  a 
priori  given  tolerance  emax.  This  defines  the  properties  of  the  ideal  mesh  that 
should  be  used. 

In  practice  we  proceed  analogously  as  in  Section  1.4.  For  each  super-element 
Sj,  there  is  a  sub-tree  in  the  quadtree  representation  of  the  current  mesh  A.  Let 

(o'  be  a  node  of  the  sub-tree  which  again  defines  a  sub-tree.  Let  now  A(co)  be 
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the  collection  of  all  elements  corresponding  to  the  terminal  nodes  of  this  sub¬ 
tree.  Then  we  can  define  the  intensity 

0)6  A(co)  u 

for  co'.  With  this  ^(co')  and  the  required  intensity  3,  the  modification  of  the  current 

2 

mesh  is  then  based  on  the  following  two  rules  for  each  2-D  domain  £2k : 

(a)  If  9? (co')  >  3  and  co'  is  a  father  node  of  an  element  co  in  the  current 
mesh,  then  the  sub-division  of  co'  is  eliminated;  that  is,  co  and  its 
brothers  are  removed  and  co’  becomes  an  element  in  a  de-refined 
mesh. 

(b)  If  SR(co')  <  3  and  co'  is  an  element  of  the  current  mesh  refined  by 
subdividing  co’  into  four  elements. 

When  a  father  node  co'  is  de-refined  and,  by  rule  (a),  becomes  an  element, 
then  rule  (a)  has  to  be  applied  again  to  its  father  node.  Similarly,  when,  by  rule 
(b),  an  element  co’  is  sub-divided  then  the  new  elements  must  also  be  checked 
by  rule  (b).  NFEARS  performs  first  the  de-refinement  with  rule  (a)  until  no  further 
mesh  modification  occurs,  then  it  proceeds  with  the  refinement  by  rule  (b),  and 
the  process  ends  when  no  more  refinement  is  needed. 
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